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COMPUTER  SIMULATION  OF  SPACE  CHARGE  EFFECTS  IN 
PARTICLE  ACCELERATORS 


I . Introduction 

Recent  studies^  have  indicated  the  feasibility  of  using  intense  en- 
ergetic heavy  ion  beams  as  the  source  of  ignition  energy  for  inertially 
confined  nuclear  fusion  reactions.  It  appears  that  heavy  ion  beams  with 
suitable  powers  and  energies  can  be  produced  using  modest  extrapolations 
of  conventional  particle  accelerator  technology.  However,  one  of  the 
major  limits  on  the  beam  intensities  which  can  be  delivered  to  a target 
IS  due  to  the  space  charge  repulsive  forces  of  the  beam  particles  them- 
selves. It  has  therefore  become  important  to  examine  the  effects  of 
these  space  charge  forces  on  particle  focussing  and  transport. 

The  transport  of  charged  particle  beams  is  well  understood  when  the 

2 

space  charge  densities  are  small.  In  this  limit,  the  particle  trajec- 
tories can  be  accurately  described  by  their  single  particle  orbits  in 
the  external  focussing  fields.  A commonly  employed  focussing  system 
uses  alternate  gradient  quadrupole  magnets.  We  will  assume  that  all 
particles  have  the  same  longitudinal  velocity  and  proceed  to  discuss  the 
transverse  behavior  of  quadrupole  systems. 

If  the  plane  perpendicular  to  the  direction  of  beam  propagation  is 
described  in  rectangular  coordinates  properly  aligned  with  respect  to 
the  quadrupole  magnet  pole  tips,  each  particle  is  subjected  to  an 

Note;  Manuscript  submitted  January  16,  1978. 
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x-directed  force  proportional  to  its  x distance  from  the  center  of  the 
system  and  a y-directed  force  proportional  to  its  y distance  from  the 
center.  If  this  force  is  focussing  in  x,  it  is  defocussing  in  y.  Each 
successive  magnet  is  then  rotated  90°  so  that  in  each  direction  a parti- 
cle is  alternatively  focussed  and  defocussed.  However,  the  net  effect 
is  focussing  in  both  directions. 


A particle  which  starts  off  center,  is  focussed  towards  the  center, 
goes  past  it,  and  eventually  returns  through  its  initial  position,  there- 
by executing  a periodic  motion.  The  fraction  of  this  period  excecuted 
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per  lens  pair,  is  a phase  shift  given  by 


Cos 


1 - n^/2 


(1) 


where  n = v^t/F,  v^t  is  the  distance  between  succeeding  lenses.  2t  is 
therefore  the  time  taken  to  traverse  a complete  lens  doublet. 

If  the  trajectory  of  a particle  in  x-p^  phase  space  is  plotted  as 
it  passes  through  successive  focussing  lenses,  in  the  non-degenerate 
case  where  i •^,  the  particle  traces  an  ellipse.  Every  other 
particle  which  has  an  initial  phase  space  location  on  this  same  ellipse 
also  remains  on  it,  so  that  an  initial  phase  space  distribution  which  is 
uniform  along  this  ellipse  remains  that  way.  A beam  which  fills  such  an 
ellipse  will  retain  a constant  cross  section  as  it  passes  through  suc- 
ceeding lens  systems  and  is  spoken  of  as  a matched  beam.  The  product 
^0^x0  ’ area  of  the  ellipse  in  phase  space  divided  by  ir, 

is  known  as  the  emittance  c. 


{ 


2 


1 


n 


Fig.  la  shows  the  x-p^  phase  space  of  such  a matched  ellipse  of 
particles  half  way  through  a focussing  lens.  Fig.  lb  shows  the  phase 
space  2/5  way  to  the  next  (defocussing)  lens.  Fig.  Ic  shows  the  phase 
space  half  way  through  the  defocussing  lens.  Fig.  Id  shows  the  phase 
space  of  the  ellipse  now  having  executed  ^/5  of  its  period  and  Fig.  le 
shows  the  ellipse  back  where  it  started.  It  should  be  noted  that  though 
the  ellipse  looks  as  it  did  when  it  started,  the  individual  particles 
have  moved  only  degrees  (in  this  case  = 90°)  around  the  ellipse. 
The  space  charge  forces  in  this  case  are  negligibly  small. 


In  the  presence  of  space  charge  forces,  the  beam  transport  system 
becomes  somewhat  more  complex.  Each  particle  orbit  is  no  longer  inde- 
pendent  of  the  others  and  the  x and  y directions  are  no  longer  decoupled. 
However,  by  assuming  a particle  distribution  function  of  the  form 


f(r,p)  = « j{x/xj2  ^ ^ (pyp^^)2  + (Py/Py^)2  - r^2j  ^ (j) 

which  is  uniformly  distributed  on  the  x-p  -y-p  hyper-ellipse, 

X y 

Kapchinskij  and  Vladimirskij^  (K-V),  have  derived  a self-consistent  set 
of  equilibria  for  focussing  systems  in  the  presence  of  space  charge.  By 
assuming  such  an  idealized  K-V  distribution  function,  the  beam  dynamics 
including  space  charge  can  be  described  by  a pair  of  coupled  ordinary 
differential  equations  for  the  two  configuration  space  axes  of  the  K-V 
hyper-ellipse.  These  may  be  written  in  the  form  of  Lambertson,  Laslett 
and  Smith^  as 
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where  9 is  a normalized  distance  along  the  beam,  u and  u are  propor- 

X y 

tional  to  the  major  axes  of  the  elliptical  beam  cross  sections  in  x and 
y.  The  focussing  term  S(9)  is  a unit  amplitude  6 function  in  the  thin 
lens  approximation.  The  dimensionless  parameter 


Q = 


4q' 


Nr 


(4) 


where  rp  is  the  classical  proton  radius  and  K=l/fL  is  the  lens  strength, 
is  a measure  of  the  strength  of  the  space  charge  effects  on  the  beam. 
Since  the  space  charge  forces  are  repulsive  and  therefore  act  against 
the  lens  focussing  forces,  the  phase  shift  per  lens  system,  u,  is  less 
than  the  value  in  the  absence  of  space  charge.  Since  it  is  a prop- 
erty of  the  K-V  distribution  that  the  space  charge  forces  are  linear,  u 
is  the  same  for  all  beam  particles  and  can  also  serve  as  a measure  of 
the  amount  of  space  charge.  It  is  significant  to  note  that  in  the  ab- 
sence of  space  charge  there  is  only  one  free  parameter.  That  is,  the 
distance  between  lenses,  which  is  related  to  the  phase  shift  per  cell 
Adding  the  space  charge  then,  the  system  of  equations  has  two  free 
parameters.  These  beam  envelope  equations  can  be  integrated  numerically 
to  obtain  the  behavior  of  this  special  class  of  space  charge  systems.^’ ^ 
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The  analysis  of  space  charge  effects  by  integrating  the  equations 
for  the  axes  of  a K-V  distribution  suffers  from  a number  of  fundamental 
limitations.  Only  uniformly  loaded  ellipses  of  the  K-V  form  are  treated. 
Therefore,  any  collective  behavior  not  characteristic  of  a uniformly 
loaded  ellipse  can  not  be  represented.  Degrees  of  freedom  expected  to 
occur  in  a real  system  such  as  deviation  from  either  uniformity  or  el- 
lipticity  can  not  be  represented.  Therefore,  even  if  the  K-V  distribu- 
tion were  a realistic  initial  distribution,  the  stability  of  this  system 
to  realistic  perturbations  can  not  be  investigated  within  the  K-V  for- 
malism. Many  physical  effects  which  can  be  of  crucial  importance  to 
successful  accelerator  design  may  be  missed  by  this  approach. 

II . The  Numerical  Model 

Standard  particle  simulation  techniques^  can  be  employed  to  inves- 
tigate finite  space  charge  beam  transport  somewhat  more  generally  than 
is  possible  by  following  the  K-V  envelope  equations.  Since  the  beam  is 
long  compared  with  its  transverse  dimensions,  non-uniformities  of  the 
beam  in  the  direction  of  propagation  are"  neglected.  By  working  in  a 
reference  frame  moving  with  the  beam,  the  evolution  of  the  beam  as  it 
propagates  in  the  z direction  translates  into  the  temporal  evolution  as 
viewed  in  the  x-y  plane.  In  this  reference  frame  the  magnet  forces  be- 
come electrostatic  forces. 

The  numerical  time  differencing  is  shown  in  Fig.  2.  The  particle 
velocities  and  posritions  are  defined  one  half  time  step  apart.  The  par- 
ticle positions  are  advanced  according  to 

Vll  = "t  * ''t+At/2 
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The  particles  at  the  new  times  are  then  accumulated  onto  a grid  to  give 
a charge  density  array.  Poisson's  equation  is  then  solved  using  fast 
Fourrier  transforms  to  give  the  electric  fields  self-consistent  with  the 
charge  density.  These  electric  fields  are  then  used  to  change  the  veloc- 
ities according  to 

''t+3At/2  ''t+&t/2  (6) 

This  numerical  algorithm  is  time  centered,  reversible  and  second  order 
accurate.  It  is  also  simple  and  fast. 


The  quadrupole  lenses  are  assumed  to  be  thin  and  implemented  numer- 
ically by  applying  to  each  particle  impulsive  accelerations,  Av  and  Av 

X y 

proportional  to  the  x and  y displacement  of  that  particle  from  the  cen- 
ter of  the  system,  v^  is  assumed  unity  for  convenience  and  the  parti- 
cles are  then  given  an  impulsive  acceleration 


Av^  = -(x-x^)6 
AVy  = (y-yo)« 


(7) 


when  passing  through  a lens  which  is  focussing  in  the  x direction  (de- 
focussing in  the  y direction).  Here  (x^,  y^)  is  the  center  of  the  sys- 
tem and  5 = 1/F,  where  F is  the  focal  length  of  the  lens.  In  the 
absence  of  space  charge  the  particles  will  then  coast  between  lenses 
according  to  Eq.  (6)  so  that  the  numerical  system  reproduces  the  thin 
lens  system  to  within  machine  roundoff. 


Because  of  the  finite  grid  on  which  the  electric  field  is  defined, 
space  charge  phenomena  can  be  resolved  only  on  the  scale  of  a grid  cell. 
Because  Poisson's  equation  is  only  solved  once  per  time  step  the  charge 


density  can  not  be  allowed  to  vary  significantly  within  a time  step. 

The  physical  implications  of  these  statements  will  be  dealt  with  in  the 
next  section. 

Ill . Tests  of  the  Numerical  Model 

Because  the  numerical  model  being  employed  follows  many  particles 
in  their  self-consistent  orbits  the  large  class  of  nonlinearities  possi- 
ble in  an  "n-body  problem,"  where  n is  several  thousand,  may  be  present. 
The  details  of  the  physical  behavior  to  be  modeled  may  therefore  be 
quite  complex.  The  numerical  system  should  do  a good  job  of  modeling 
the  "important"  behavior,  where  important  is  defined  in  terms  of  the  be- 
havior of  the  average  system  parameters.  For  example,  if  the  physical 
system  is  unstable  to  a perturbation  whose  characteristic  wavelength  is 
smaller  than  a cell  size,  this  behavior  may  be  missed  entirely,  or  at 
least  poorly  represented,  because  of  the  inability  of  the  numerical  sys- 
tem to  resolve  phenomena  on  this  scale.  The  same  may  be  true  for  per- 
turbations with  growth  times  much  smaller  than  a timestep.  Lacking  a 
complete  theory  of  the  system  behavior,  such  questions  must  be  resolved 
experimentally  by  varying  the  numerical  parameters  until  a regime  is 
found  where  the  system  behavior  is  insensitive  to  numerical  parameters. 
From  similar  simulations  in  plasma  physics  where  particle  dynamics  are 
dominated  by  the  collective  nature  of  the  electric  fields,  it  may  be 
expected  that  the  time  step  must  be  of  order  the  plasma  oscillation  time 
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Preliminary  designs  for  transport  systems^  have  assumed  a rough 
balance  between  magnetic  focussing  forces  and  space  charge  defocussing 
forces.  Accordingly,  numerical  tests  have  been  run  on  a transport 
system  which  is  in  this  range  of  parameters. 


Magnet  strength  in  the  test  case  has  been  chosen  to  correspond  to  a 
phase  shift  per  doublet  of  90°  as  defined  in  Eq.  1.  In  the  absence  of 
space  charge,  as  in  the  case  which  was  shown  in  Fig.  1,  it  is  possible 
to  calculate  the  initial  conditions  for  a matched  system  analytically. 
However,  in  the  presence  of  space  charge,  it  is  most  convenient  to  cal- 
culate the  initial  condition  for  a matched  system  iteratively.  Succes- 
sive guesses  are  obtained  by  following  the  system  through  a doublet  lens 
system,  from  t=o  to  t=2T.  The  x extent  of  the  next  guess  is  chosen  as 
the  average  of  x extent  initially  and  after  one  doublet  holding  the  emit- 
tance,  x p^,  constant.  Similarly,  the  y extent  at  t=2T  is  used  to  cal- 
culate a new  averaged  y dimension.  This  procedure  usually  converges  in 
just  a few  iterations.  Fig.  3 shows  the  initial  conditions  thus  found 
at  the  center  of  a lens  which  is  focussing  in  the  x direction,  for  a 
slightly  mismatched  system  with  a phase  shift  per  doublet  of  = ir/2. 

A mismatch  of  about  10?o  was  chosen  so  the  system  would  evolve  fairly 
rapidly  and  thereby  test  the  numerics.  Four  projections  of  the  density 
in  x-p^-y-Py  space  are  shown  for  2000  sample  particles.  The  distances 
labelled  on  the  axes  are  in  terms  of  program  grid  points  and  velocities 
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are  in  units  of  cbIIs/t.  The  initial  space  charge  density  is  such  that 
the  Q defined  in  Eq.  h is  equal  to  3.9.  In  this  case  the  phase  shift 
per  doublet  is  about  30°.  The  particle  positions  and  momenta  are  set  up 
using  the  formulae 

X = X cos  hr-) 

0 Zir 

Px  = Pxo 

y = Yq  cos  (-|^) 

p = p /T^  sin  (-1—) 

"^yo  ^ 2n 

where  q,  r,  and  s are  uniformly  distributed  random  numbers  between  zero  ' 

and  unity. 

The  Poisson  solver  uses  Fourier  transforms  and  periodic  boundary 
conditions  to  solve  for  the  electric  field  from  the  particle  charge  dis- 
tribution. Since  it  is  impractical  to  follow  the  trajectories  of  as 
many  particles  as  are  actually  found  in  a beam,  anomolously  high  fluc- 
tuations in  density  occur  in  a particle  code  simulation.  In  addition, 
following  only  a finite  number  of  modes  leads  to  a band-limited  situa- 
tion in  which  short  wavelength  overshoots  can  occur  in  the  vicinity  of 
sharp  changes  in  charge  density.  Both  these  effects  can  be  mitigated 
somewhat  by  considering  "macro-particles"  with  a Gaussian  charge  den-  j 

sity.^  When  this  is  done,  the  points  plotted  as  in  Fig.  3 are  really 
the  centers  of  these  "macroparticle"  charge  clouds  rather  than  point 
particles. 
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Fig.  4a  is  a plot  of  vs  x for  a y coordinate  1/2  cell  off  the 
center  of  a 64  x 64  system  whose  initial  conditions  are  similar  to  those 
shown  in  Fig.  3.  Fig.  4b  is  a plot  of  the  same  fields  with  Gaussian 
smoothing  applied  to  the  system  so  that  the  particle  "size"  between  1/e 
points  on  the  Gaussian  is  one  cell.  Note  that  is  really  the  average 
field  applied  to  the  smoothed  cloud,  i.e.,  the  acceleration  field.  The 
smoothed  fields  are  closer  to  linear  near  the  center  but  the  peaks  at 
the  edge  of  the  beam  are  somewhat  rounded.  Figs.  5 and  6 are  a compari- 
son of  the  phase  space  of  two  systems  after  traversing  16  focussing  mag- 
net pairs  with  and  without  smoothing.  The  comparison  time  was  chosen 
because  of  the  fine  structure  at  this  time  which  ought  to  accentuate 
numerical  differences.  Figs.  5 and  6 are  phase  space  plots  which  corre- 
spond to  the  same  smoothing  as  in  Figs.  4a  and  4b  respectively.  Though 
there  are  minor  differences,  they  do  not  appear  to  be  significant  and 
they  are  small  compared  with  differences  which  arise  from  changing  the 
number  of  grid  cells  and  which  will  be  dealt  with  later.  In  both  Figs. 

5 and  6,  2048  particles  are  shown  out  of  the  4096  used  in  the  run. 

Other  parameters  are  At  = .2t,  64  x 64  cells  and  (UpX  = 1.2.  Since  the 
smoothing  does  not  have  a major  effect  on  the  short  time  system  evolu- 
tion and  does  not  seem  to  negatively  affect  evolution  of  the  beam  edge, 
there  is  a net  advantage  to  employing  it  for  long  times  in  order  to 
avoid  anomolous  fluctuation  effects. 

In  simulating  a collisionless  plasma,  it  is  usually  impractical  to 
represent  anything  like  the  number  of  particles  in  a real  system.  Con- 
sequently the  number  of  particles  per  Debye  volume  is  low.  The  simula- 
tion may  then  represent  the  collisionless  behavior  only  for  a number  of 
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plasma  periods  proportional  to  the  reciprocal  of  the  actual  number  of 
simulation  particles  in  a Debye  volume.  Fig.  7 is  a phase  space  plot 
from  a calculation  with  conditions  as  in  Fig.  6 but  with  double  the  num- 
ber (8192)  of  simulation  particles.  While  the  fine  structure  changes 
slightly,  mainly  in  the  y“Py  view,  the  computed  results  are  essentially 
independent  of  the  number  of  simulation  particles;  the  small  residual 
differences  are  due  as  much  to  different  statistics  in  the  initial  dis- 
tribution as  to  fluctuations  in  the  fields  during  the  run.  The  differ- 
ence is  made  more  noticeable  in  Fig.  8 by  again  doubling  the  number  of 
particles  to  16384.  Here,  too,  however,  it  is  felt  that  any  differences 
are  due  to  different  statistics. 

i 

I 

I The  Fourier  Transform  Poisson  solver  employed  uses  periodic  bound- 

ary conditions.  This  means  that  the  neutralizing  charges,  instead  of 
being  at  infinity  or  in  finite  images  are,  uniformly  distributed  through- 
out the  system.  Furthermore  each  beam  particle  has  an  image  of  the  same 
sign  and  size  one  system  length  away  in  each  direction.  Such  a computa- 
tional system  clearly  reproduces  the  infinite  system  when  the  boundary 
is  far  enough  away. 

Fig.  9 shows  results  from  a calculation  using  a 128  by  128  finite 
difference  mesh  where  the  beam  occupies  only  one  quarter  of  the  system 
area  of  the  previous  runs.  That  is,  measured  in  beam  size,  the  bounda- 
ries are  twice  as  far  out  as  in  the  previous  runs.  The  numerics  are 
otherwise  identical  to  Fig.  6.  The  behavior  of  the  two  systems  is  again 
quite  similar,  showing  that  the  boundary  conditions  are  not  having  a 
major  effect. 
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We  have  thus  far  seen  that,  when  run  through  16  magnet  pairs,  the  system 
behavior  is  relatively  insensitive  to  varying  the  particle  smoothing, 
the  number  of  particles  and  the  distance  to  the  boundaries. 

A comparison  of  Fig.  6 with  Fig.  10  shows  the  differences  in  be- 
havior caused  by  changes  in  the  resolution.  Fig.  6 uses  twice  as  many 
cells  in  the  description  of  the  electric  fields  of  the  beam  and  fine 
structure  has  appeared  which  was  not  present  in  the  coarser  grid  run  of 
Fig.  10.  A comparison  of  Figs.  9 and  11  shows  the  result  of  doubling 
the  resolution  again.  Fig.  9 is  a 64  x 64  system  and  Fig.  11  is  128  x 
128.  Each  run  has  16K  particles.  (The  number  of  sample  particles  for 
plotting  is  still  2K.)  Though  the  difference  is  not  as  great  as  the 
original  doubling,  more  fine  structure  is  evident. 

Finally,  Fig.  12  shows  the  effect  of  taking  twice  as  many  timesteps 
(20)  per  magnet  period.  The  parameters  are  otherwise  the  same  as  the 
run  in  Fig.  11.  A comparison  of  the  two  pictures  shows  that  the  two  are 
extremely  close. 

It  is  evident  from  these  series  of  tests  that  the  numerics  are  in- 
deed capable  of  affecting  the  fine  structure  of  the  phase  space  distri- 
bution of  the  beam  particles.  The  important  question  is  whether  differ- 
ences in  the  fine  structure  have  any  important  physical  significance. 
Though  the  question  can  only  be  answered  definitively  when  the  physics 
of  space  charge  effects  in  particle  transport  systems  is  well  understood, 
there  are  several  reasons  for  confidence  that  the  numerical  results  will 
be  reliable. 
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The  fine  structure  observed  in  these  runs  is  probably  a consequence 
of  the  sharp  edge  of  the  distribution  function.  This  particular  initial 
distribution  was  chosen  because  the  K-V  distribution  is  well  understood 
theoretically.  It  is  not  a physically  realistic  distribution  because  of 
the  very  sharp  edges.  A more  physical  distribution  would  have  diffuse 
edges  and  the  fine  structured  behavior  would  then  be  washed  out.  Some 
evidence  has  already  been  obtained  that  this  is  the  case.  In  fact,  it 
appears  as  if  a steady  state  can  be  reached  which  is  independent  of  the 
numerical  parameters  and  the  details  of  the  intermediate  turbulent 
structure. 

That  the  numerical  system  can  do  as  good  a job  as  it  does  in  exam- 
ining the  behavior  of  the  sharp  edged  K-V  distribution  is  encouraging 
since  this  type  of  system  is  generally  the  most  difficult  one  to  treat 
numerically.  It  seems  clear  that  the  similarities  in  behavior  over  a 
large  range  of  numerical  parameters  indicate  that  reliable  data  on 
focussed  beam  systems  can  be  inferred  from  this  approach.  Furthermore, 
because  of  the  flexibility  of  this  numerical  approach  it  is  easily  pos- 
sible to  include  many  complexities  such  as  finite  length  lenses  and  lens 
non-linearities  which  can  give  valuable  insight  into  design  of  practical 
transport  systems.  Initial  simulations  to  gain  insight  into  the  physics 
of  beam  transport  systems  in  the  presence  of  space  charge  should  concen- 
trate on  simple  systems  and  these  are  currently  underway. 
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Fig.  1 — x-Pjf  phase  space  plot  showing  the  evolution  of  a matched  phase 
space  ellipse  of  particles  as  it  passes  through  a lens  system  with  a phase 
shift  of  90  per  doublet.  Fig.  la  shows  the  initial  distribution  halfway 
through  a focussing  lens.  Fig.  lb  shows  the  distribution  2/5  of  the  dis- 
tance to  the  next  lens.  Fig.  Ic  shows  the  distribution  halfway  throu^  a 
defocussing  lens.  Fig.  Id  is  4/5  throu^  the  doublet.  Finally,  Fig.  le  is 
one  complete  lens  period  away  showing  a return  to  the  initial  elliptical 
shape. 
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Fig.  2 — Time  centering  of  the  numerical  system  showing  how  the 
system  quantities  are  advanced  through  a time  step 


Fig.  3 — Initial  phase  space  ellipses  of  a slightly  mismatched  Kapchinskij- 
Vladimirkij  distribution  with  a magnet  tune  of  90°  per  doublet  and  e- 
nough  space  charge  to  bring  the  tune  down  to  30°.  Units  of  length  are 
in  cells  and  units  of  velocity  are  in  cells/r  where  r is  the  time  between  suc- 
cessive magnets.  Center  of  the  beam  is  the  center  of  the  numerical  system 
so  that  the  system  size  (128  x 128)  can  be  inferred  from  the  labels  on  the  . 
configuration  space  axes.  2048  sample  particles  are  plotted. 
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Fig.  5 — Phase  space  of  a 64  x 64  system  showing  2048  sample  particles 
of  4096  in  the  system  after  16  magnet  pairs  of  160  time  steps,  (r  = 32) 
No  smoothing  was  used. 
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Fig.  6 — Phase  space  of  a 64  x 64  system  with  a 1/e  width  of  one  cell 
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Fig.  8 — Phase  space  of  a 64  x 64  system  with  the  number  of  particles 
again  doubled  (to  16284) 
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Fig.  9 — Phase  space  of  a 128  x 128  system  with  4096  particles. 
The  resolution  describing  the  beam  is  unchanged  but  the  bounda- 
ries are  twice  as  far  away. 
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Fig.  11  — Phase  space  of  2048  sample  particles  on  a system  with 
sixteen  thousand  particles  and  128  x 128  grid  points  showing  the 
effect  of  the  greater  resolution  on  the  fine  structure  of  the  system 
evolution 


